/*
 * Float Array Validation with Error Tolerance
 * Implements comprehensive validation framework for GPU correctness verification
 */

#include <cuda_runtime.h>
#include <math.h>
#include <cstdio.h>
#include <cstdlib.h>
#include <string>
#include <vector>

// Validation configuration structure
struct ValidationConfig {
    float tolerance = 1e-6f;
    bool verbose = false;
    bool check_nan = true;
    bool check_infinity = true;
    bool check_denormals = true;
    int max_error_report = 10;
    bool generate_summary = true;
};

// Validation result structure
struct ValidationResult {
    bool passed = false;
    int total_elements = 0;
    int error_count = 0;
    float max_error = 0.0f;
    int max_error_index = -1;
    int nan_count = 0;
    int inf_count = 0;
    int denormal_count = 0;
    std::vector<int> error_indices;
    std::vector<float> error_values;
    std::vector<float> expected_values;
    
    void add_error(int index, float computed, float expected) {
        error_indices.push_back(index);
        error_values.push_back(computed);
        expected_values.push_back(expected);
        
        float diff = fabsf(computed - expected);
        if (diff > max_error) {
            max_error = diff;
            max_error_index = index;
        }
    }
    
    void print_summary() const {
        printf("Validation Summary:\n");
        printf("  Total elements: %d\n", total_elements);
        printf("  Errors: %d\n", error_count);
        printf("  Pass rate: %.2f%%\n", 
               (total_elements > 0) ? (100.0f * (total_elements - error_count) / total_elements) : 0.0f);
        
        if (max_error_index >= 0) {
            printf("  Max error: %.6f at index %d\n", max_error, max_error_index);
        }
        
        if (nan_count > 0) printf("  NaN values: %d\n", nan_count);
        if (inf_count > 0) printf("  Infinity values: %d\n", inf_count);
        if (denormal_count > 0) printf("  Denormal values: %d\n", denormal_count);
    }
};

// GPU kernel for error checking
__global__ void gpu_error_checking_kernel(
    const float* __restrict__ computed,
    const float* __restrict__ expected,
    int* __restrict__ error_flags,
    float* __restrict__ error_values,
    int* __restrict__ error_indices,
    int* __restrict__ error_count,
    int n,
    float tolerance
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (idx < n) {
        float diff = fabsf(computed[idx] - expected[idx]);
        
        // Check for NaN and infinity
        bool is_nan_error = isnan(computed[idx]) != isnan(expected[idx]);
        bool is_inf_error = isinf(computed[idx]) != isinf(expected[idx]);
        
        if (is_nan_error || is_inf_error || diff > tolerance) {
            int error_idx = atomicAdd(error_count, 1);
            error_flags[idx] = 1;
            
            if (error_idx < 1000) {  // Limit storage
                error_values[error_idx * 2] = computed[idx];
                error_values[error_idx * 2 + 1] = expected[idx];
                error_indices[error_idx] = idx;
            }
        }
    }
}

// CPU validation class
class FloatArrayValidator {
private:
    ValidationConfig config;
    
public:
    FloatArrayValidator(const ValidationConfig& cfg = ValidationConfig()) : config(cfg) {}
    
    // Main validation function
    ValidationResult validate(const float* computed, const float* expected, int size) {
        ValidationResult result;
        result.total_elements = size;
        result.error_count = 0;
        
        for (int i = 0; i < size; i++) {
            bool has_error = false;
            
            // Check for NaN
            if (config.check_nan) {
                bool computed_nan = isnan(computed[i]);
                bool expected_nan = isnan(expected[i]);
                
                if (computed_nan != expected_nan) {
                    result.nan_count++;
                    has_error = true;
                }
            }
            
            // Check for infinity
            if (config.check_infinity) {
                bool computed_inf = isinf(computed[i]);
                bool expected_inf = isinf(expected[i]);
                
                if (computed_inf != expected_inf) {
                    result.inf_count++;
                    has_error = true;
                }
            }
            
            // Check for denormals
            if (config.check_denormals) {
                bool computed_denorm = (fabsf(computed[i]) < 1e-38f && computed[i] != 0.0f);
                bool expected_denorm = (fabsf(expected[i]) < 1e-38f && expected[i] != 0.0f);
                
                if (computed_denorm != expected_denorm) {
                    result.denormal_count++;
                    has_error = true;
                }
            }
            
            // Check numerical error
            if (!isnan(computed[i]) && !isnan(expected[i]) &&
                !isinf(computed[i]) && !isinf(expected[i])) {
                
                float diff = fabsf(computed[i] - expected[i]);
                
                // Handle relative error for large values
                float relative_error = 0.0f;
                if (fabsf(expected[i]) > 1e-10f) {
                    relative_error = diff / fabsf(expected[i]);
                } else {
                    relative_error = diff;  // Absolute error for small values
                }
                
                if (diff > config.tolerance && relative_error > config.tolerance) {
                    result.add_error(i, computed[i], expected[i]);
                    result.error_count++;
                    has_error = true;
                }
            }
            
            if (has_error) {
                if (config.verbose && result.error_count <= config.max_error_report) {
                    printf("Error at index %d: computed=%.8f, expected=%.8f\n", 
                           i, computed[i], expected[i]);
                }
            }
        }
        
        result.passed = (result.error_count == 0);
        
        if (config.generate_summary) {
            result.print_summary();
        }
        
        return result;
    }
    
    // GPU-accelerated validation
    ValidationResult validate_gpu(const float* d_computed, const float* d_expected, int size) {
        ValidationResult result;
        result.total_elements = size;
        
        // Allocate GPU memory for error checking
        int* d_error_flags = nullptr;
        int* d_error_count = nullptr;
        float* d_error_values = nullptr;
        int* d_error_indices = nullptr;
        
        cudaMalloc(&d_error_flags, size * sizeof(int));
        cudaMalloc(&d_error_count, sizeof(int));
        cudaMalloc(&d_error_values, 2000 * sizeof(float));  // Max 1000 errors
        cudaMalloc(&d_error_indices, 1000 * sizeof(int));
        
        cudaMemset(d_error_flags, 0, size * sizeof(int));
        cudaMemset(d_error_count, 0, sizeof(int));
        
        // Launch kernel
        int threads_per_block = 256;
        int blocks = (size + threads_per_block - 1) / threads_per_block;
        
        gpu_error_checking_kernel<<<blocks, threads_per_block>>>(
            d_computed, d_expected, d_error_flags, d_error_values,
            d_error_indices, d_error_count, size, config.tolerance
        );
        
        cudaDeviceSynchronize();
        
        // Copy results back
        int h_error_count = 0;
        cudaMemcpy(&h_error_count, d_error_count, sizeof(int), cudaMemcpyDeviceToHost);
        
        result.error_count = h_error_count;
        
        if (h_error_count > 0) {
            // Copy error details
            std::vector<float> h_error_values(h_error_count * 2);
            std::vector<int> h_error_indices(h_error_count);
            
            cudaMemcpy(h_error_values.data(), d_error_values, 
                      h_error_count * 2 * sizeof(float), cudaMemcpyDeviceToHost);
            cudaMemcpy(h_error_indices.data(), d_error_indices, 
                      h_error_count * sizeof(int), cudaMemcpyDeviceToHost);
            
            for (int i = 0; i < h_error_count; i++) {
                result.add_error(
                    h_error_indices[i],
                    h_error_values[i * 2],
                    h_error_values[i * 2 + 1]
                );
            }
        }
        
        result.passed = (result.error_count == 0);
        
        // Cleanup
        cudaFree(d_error_flags);
        cudaFree(d_error_count);
        cudaFree(d_error_values);
        cudaFree(d_error_indices);
        
        if (config.generate_summary) {
            result.print_summary();
        }
        
        return result;
    }
    
    // Batch validation for multiple arrays
    bool validate_batch(const float** computed_arrays, const float** expected_arrays,
                       int array_size, int num_arrays, ValidationResult* results = nullptr) {
        bool all_passed = true;
        
        for (int i = 0; i < num_arrays; i++) {
            ValidationResult result = validate(computed_arrays[i], expected_arrays[i], array_size);
            
            if (results) {
                results[i] = result;
            }
            
            if (!result.passed) {
                all_passed = false;
                printf("Array %d validation failed (%d errors)\n", i, result.error_count);
            }
        }
        
        return all_passed;
    }
    
    // Statistical analysis
    void calculate_statistics(const float* computed, const float* expected, int size,
                            float* mean_error, float* max_error, float* stddev_error) {
        double sum_error = 0.0;
        double sum_squared_error = 0.0;
        double max_err = 0.0;
        
        for (int i = 0; i < size; i++) {
            if (!isnan(computed[i]) && !isnan(expected[i]) &&
                !isinf(computed[i]) && !isinf(expected[i])) {
                
                float diff = fabsf(computed[i] - expected[i]);
                sum_error += diff;
                sum_squared_error += diff * diff;
                max_err = fmax(max_err, diff);
            }
        }
        
        *mean_error = (float)(sum_error / size);
        *max_error = (float)max_err;
        *stddev_error = (float)sqrt(sum_squared_error / size - (*mean_error) * (*mean_error));
    }
    
    // Configuration setters
    void set_tolerance(float tolerance) { config.tolerance = tolerance; }
    void set_verbose(bool verbose) { config.verbose = verbose; }
    void set_max_error_report(int max_report) { config.max_error_report = max_report; }
};

// C-style interface functions
extern "C" {
    // Simple validation
    int validate_float_arrays(const float* computed, const float* expected, int size,
                            float tolerance, int* error_count) {
        ValidationConfig config;
        config.tolerance = tolerance;
        config.verbose = false;
        
        FloatArrayValidator validator(config);
        ValidationResult result = validator.validate(computed, expected, size);
        
        if (error_count) *error_count = result.error_count;
        return result.passed ? 1 : 0;
    }
    
    // Detailed validation
    void validate_arrays_detailed(const float* computed, const float* expected, int size,
                                float tolerance, ValidationResult* result) {
        ValidationConfig config;
        config.tolerance = tolerance;
        config.verbose = true;
        
        FloatArrayValidator validator(config);
        *result = validator.validate(computed, expected, size);
    }
    
    // Batch validation
    int validate_batch_arrays(const float** computed, const float** expected,
                            int array_size, int num_arrays, float tolerance) {
        ValidationConfig config;
        config.tolerance = tolerance;
        config.verbose = false;
        
        FloatArrayValidator validator(config);
        return validator.validate_batch(computed, expected, array_size, num_arrays) ? 1 : 0;
    }
    
    // Statistics calculation
    void calculate_validation_stats(const float* computed, const float* expected, int size,
                                  float* mean_error, float* max_error, float* stddev_error) {
        ValidationConfig config;
        FloatArrayValidator validator(config);
        validator.calculate_statistics(computed, expected, size, mean_error, max_error, stddev_error);
    }
}